/***************************************************************
                       score2.c

only routines necessary for parallel
**************************************************************/
#include <stdio.h>
#include <malloc.h>
#include <float.h>
#include <math.h>

#include <assert.h>
#include <netinet/in.h>
#include "serial.h"

#include "clam.h"

#define FAST_SUL_ITERATIONS 3

#define graphQuery messQuery

void CpMdisplay(), ZautoCpM();
int  ZscoreCpM();
BOOL ZboolCpM();

static CLONE my_clone;
extern void scanPzA(), scanPzB(), scanPz2(), EvalCtgDisplay(), CBCtgDisplay(), ClamCtgDisplay2();
extern void ZclearHigh(), ZhighCin();
extern int Zget_tol();
extern void ClamShare();
extern int NoPop;
extern int validDistribution;
extern float *distrib;

/* sness */
static void fast_sulston_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result);
static void pure_sulston_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result);

void Zsulston_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result);

void Zgetmatch_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result);
int ZscoreCpM_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result);
int Zcnt_shared_markers_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result);
void Zmott_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result);

extern int get_precomp_value(int b1, int b2, int m, double* result);

BOOL ZboolCpM();
int Zcnt_shared_markers();
void Read_unc(), set_Zunc(), ZbuildCpM();

/*************************************************************
                     CONTIGC
Finds the probability of coincidence
C(n, k) = (N over K) = bionomial coeffient = N!/K!(N-K)!

sum of N to nbandL C(nbandL, N) * (1 - p)^N * p^(nbandL-N)
where p = (1 - (2*tol/gel-len)) ^ nbandH

**************************************************************/
#define abs(x) ((x) >= 0 ? (x) : -(x))

/* */
/* CzFastCache stores coordinate and band information */
/* */

/* Cache data for loadCzfast */
static unsigned short*  coords;
static int     coordsSize;
static int*    coordOffset;
static unsigned short*  bandCount;
static int     nclones = -1;
static int     minVal = 0;
static int     maxVal = 0;

int releaseCzFastCache()
{
	if (nclones != -1)
	{
		free(coords);
		free(coordOffset);
		free(bandCount);
	}
	
	nclones = -1;
	minVal = 0;
	maxVal = 0;

	return 0;
}

int writeLoadCzFastCache(int socket)
{

     
	union chunk out[4];
	int32_t* buf;
	int i;
	
	assert(sizeof(int32_t)     == 4);
	assert(sizeof(union chunk) == 4);
	assert(sizeof(out)         == 4*4);
	
	out[0].i = htonl(coordsSize);
	out[1].i = htonl(nclones);
	out[2].i = htonl(minVal);
	out[3].i = htonl(maxVal);
	
	printf("\nclone cache coordsSize:%d nclones:%d  minVal:%d maxval:%d\n",coordsSize,nclones,minVal,maxVal);
	

	if (write(socket, &out, sizeof(out)) != sizeof(out)) return -1;
	
	//sleep(10);

	if (!(buf = (int32_t*)malloc(4 * coordsSize)))
	{
	  printf("out of memory writing cache to tcp stream.\n");fflush(stdout);
		perror("out of memory writing cache to tcp stream.");
		exit(-1);
	}
	
	for (i = 0; i < coordsSize; i++)
		buf[i] = htonl(coords[i]);
	
	if (write(socket, buf,  4*coordsSize) != 4*coordsSize)
	{
		free(buf);
		return -1;
	}
	
	free(buf);
	
	if (!(buf = (int32_t*)malloc(4 * nclones)))
	{
	  printf("out of memory writing cache to tcp stream.\n");fflush(stdout);
		perror("out of memory writing cache to tcp stream.");
		exit(-1);
	}
	
	for (i = 0; i < nclones; i++)
		buf[i] = htonl(coordOffset[i]);
	
	if (write(socket, buf, 4*nclones) != 4*nclones)
	{
		free(buf);
		return -1;
	}
	
	for (i = 0; i < nclones; i++)
		buf[i] = htonl(bandCount[i]);
	
	if (write(socket, buf, 4*nclones) != 4*nclones)
	{
		free(buf);
		return -1;
	}
	
	//sleep(10);
	free(buf);
	
	return 0;
}

int readLoadCzFastCache(int socket)
{
	union chunk in[4];
	int i;
	int32_t* buf;

	

	assert(sizeof(int32_t)     == 4);
	assert(sizeof(union chunk) == 4);
	assert(sizeof(in)          == 4*4);


	if (readLoop(socket, &in, sizeof(in)) != sizeof(in))
		return -1;

	if (nclones != -1)
	{
		free(coords);
		free(coordOffset);
		free(bandCount);
	}
	
	coordsSize = ntohl(in[0].i);
	nclones    = ntohl(in[1].i);
	minVal     = ntohl(in[2].i);
	maxVal     = ntohl(in[3].i);


	/* printf("Received clones :%d coordsize: %d\n", nclones,coordsSize); fflush(stdout); */
	
	coords      = (unsigned short*)malloc(sizeof(unsigned short) * coordsSize);
	coordOffset = (int*)  malloc(sizeof(int)   * nclones);
	bandCount   = (unsigned short*)malloc(sizeof(unsigned short) * nclones);
	
	buf = (int32_t*)malloc(4 * MaX(coordsSize, nclones));
	
	if (!coords || !coordOffset || !bandCount || !buf)
	{
	  printf("\nOut of memory loading cache from TCP.\n\n");fflush(stdout);
		perror("Out of memory loading cache from TCP.\n");
		exit(-1);
	}
	
	readLoop(socket, buf, coordsSize * 4);
	for (i = 0; i < coordsSize; i++)
		coords[i] = ntohl(buf[i]);
	
	readLoop(socket, buf, nclones * 4);
	for (i = 0; i < nclones; i++)
		coordOffset[i] = ntohl(buf[i]);
	
	readLoop(socket, buf, nclones * 4);
	for (i = 0; i < nclones; i++)
		bandCount[i] = ntohl(buf[i]);
	
	free(buf);

	return 0;
}

/* This builds the cache */
int rebuildLoadCzFastCache()
{
	int clone;
	int band;
	int bandEnd;
	int writer;
	
	
	unsigned short bandData;
	
	CLONE* clp;

	//printf("\n\n##################Building cache...\n"); fflush(stdout);
	

	
/* 	printf("\n"); */
/* 	printf("Building cache table for loadCzfast:\n"); */
	
	if (nclones != -1)
	{
		printf("  - Freeing old loadCzfast cache data.\n");
		free(coords);
		free(coordOffset);
		free(bandCount);
	}
	
	/* Update the globals */
	nclones = arrayMax(acedata);
	minVal  = Pz.minVal;
	maxVal  = Pz.maxVal;
	
	/* Determine how big the coordinate cache needs to be */
/* 	printf("  - Counting coordinates: "); */
	fflush(stdout);
	
	coordsSize = 0;
	for (clone = 0; clone < nclones; clone++)
	{
		clp = arrp(acedata, clone, CLONE);
		
		if (clp->fp == NULL) continue;
		if (clp->fp->b2 == 0) continue;
		
		band    = clp->fp->b1 - 1; /* not sure why we need this -1 */
		bandEnd = clp->fp->b2 + band;
		
		for (; band < bandEnd; band++)
		{
			bandData = arr(bands, band, unsigned short);
			if (bandData > minVal && bandData < maxVal)
				coordsSize++;
		}
	}
	
	printf("%i\n", coordsSize);
	
	/* We can now allocate data */
	printf("  - Allocating cache:");
	fflush(stdout);
	
	coords      = (unsigned short*)malloc(sizeof(unsigned short) * coordsSize);
	printf(" coords=%i", sizeof(unsigned short) * coordsSize);
	fflush(stdout);
	
	coordOffset = (int  *)malloc(sizeof(int  ) * nclones);
	printf(" coordOffset=%i", sizeof(int) * nclones);
	fflush(stdout);
	
	bandCount   = (unsigned short*)malloc(sizeof(unsigned short) * nclones);
	printf(" bandCount=%i", sizeof(unsigned short) * nclones);
	fflush(stdout);

	printf(".\n");
	
	if (!coords || !coordOffset || !bandCount)
	{
		printf("OUT OF MEMORY!\n");
		exit(-1);
	}
	
	/* Fill in the cache */
	printf("  - Filling the cache: ");
	fflush(stdout);
	writer = 0;
	for (clone = 0; clone < nclones; clone++)
	{
		clp = arrp(acedata, clone, CLONE);
		
		if (clp->fp == NULL || clp->fp->b2 == 0)
		{
			bandCount  [clone] = 0;
			coordOffset[clone] = -1;
			continue;
		}
		
		band    = clp->fp->b1 - 1;
		bandEnd = clp->fp->b2 + band;
		
		coordOffset[clone] = writer;
		
		for (; band < bandEnd; band++)
		{
			bandData = arr(bands, band, unsigned short);

			if (bandData > minVal && bandData < maxVal)
			{
				coords[writer++] = bandData;
			}
		}
		
		bandCount[clone] = writer - coordOffset[clone];
	}
	printf("done\n");
	
	printf("  - Cache built successfully.\n");
}


/*********************************************************
                    DEF: loadCzfast
**********************************************************/
int loadCzfast(struct tmpCz *Cx, int cin)
{
	if (arrayMax(acedata) != nclones || minVal != Pz.minVal || maxVal != Pz.maxVal)
		rebuildLoadCzFastCache();
		  
	
/* 	result->match = result->olap = 0; */
/* 	result->prob = 0.0; */
	
	if (coordOffset[cin] == -1)
	{
		return 0;
	}
	
	Cx->coords = coords + coordOffset[cin];
	Cx->nbands = bandCount[cin];
	Cx->cin    = cin;
	
	return 1;
}


/*********************************************************
                    DEF: copyLocalCz
**********************************************************/
void copyLocalCz(struct tmpCz* local)
{
	int i;
	
	for (i = 0; i < local->nbands; i++)
	{
		local->localCoords[i] = local->coords[i];
	}
	
	local->coords = &local->localCoords[0];
}


/*********************************************************8
                    DEF: loadCz
**********************************************************/
int loadCz(Cx, cin)
struct tmpCz *Cx;
int cin;
{
register int i, off, b, j;

  Sz.msg1[0] = Sz.msg2[0] = '\0';
  Sz.match = Sz.olap = 0;
  Sz.prob = 0.0;

  my_clone = arr(acedata,cin,CLONE);
  if (my_clone.fp == NULL) return 0;
  if (my_clone.fp->b2 == 0) return 0; /* FIX 28may99 */

  if (my_clone.mattype==PARENT) strcpy(Cx->parent," canon"); 
  else if (my_clone.mattype==PSPARENT) strcpy(Cx->parent," canon+"); 
  else strcpy(Cx->parent, my_clone.match);

  if (my_clone.mattype==EXACT) strcpy(Cx->btype,"exact "); 
  else if (my_clone.mattype==APPROX) strcpy(Cx->btype,"approx"); 
  else if (my_clone.mattype==PSEUDO) strcpy(Cx->btype,"pseudo"); 
  else strcpy(Cx->btype,"      "); 

  strcpy(Cx->clone,my_clone.clone);
  strcpy(Cx->fpnum,my_clone.fp->fpchar);

  Cx->left = my_clone.x;
  Cx->right = my_clone.y;
  Cx->len = ( Cx->right - Cx->left);
  Cx->ctg = my_clone.ctg;

  if (merging) {
     if (Cx->ctg == mergecontig2) {
         off = (contigs[mergecontig1].right - 
                contigs[mergecontig2].left) + startpt;
         Cx->left += off;
         Cx->right += off;
     }
  }

  Cx->start = my_clone.fp->b1;

  /* Removed from wterpstra's code, probably with good reason due to loadCzfast() below */
  Cx->cin = cin;



  if (my_clone.fp->b2 > NBANDS) printf("Bad number of bands %s %d\n",Cx->clone, 
         my_clone.fp->b2);

  /* Removed from wterpstra's code, probably with good reason due to loadCzfast() below */
  
  return loadCzfast(Cx, cin);


}
/*********************************************************8
                    DEF: loadSizeCz
**********************************************************/
int loadSizeCz(Cx, cin)
struct tmpCz *Cx;
{
CLONE *clone;
int tmp[NBANDS];
FILE *bandsFilePtr;
char bandsFileName[1024];
char name[CLONE_SZ];
char line[100];
int  temp;
int i,  nbands;

   if (Proj.variable) return 1; /* already have size */
   if (Pz.useSz<=0) return 0; 

   clone = arrp(acedata,cin,CLONE);
   if (strlen(dirName) > 0)
       sprintf(bandsFileName, "%s/Sizes/%s.sizes", dirName, clone->fp->gelname);
   else  sprintf(bandsFileName, "Sizes/%s.sizes",  clone->fp->gelname);

   if ((bandsFilePtr = fopen(bandsFileName, "r"))==NULL)  {
       printf("Warning: could not find %s\n",bandsFileName);
       return 0;
   }
   nbands=0;
   while (1) {
      if (!fgets(line, 80, bandsFilePtr)) {
         printf("Warning: could not find %s in %s\n",clone->clone, bandsFileName);
         return 0;
      }
      sscanf(line, "%s", name);
      if (strcmp(name,clone->clone)!=0) continue;
      fscanf(bandsFilePtr, "%i", &temp);
      while (temp != -1) {
          tmp[nbands++] = temp;
          fscanf(bandsFilePtr, "%i", &temp);
      }
      tmp[nbands] = -1;
      break;
    }
    
  
    Cx->coords = &Cx->localCoords[0]; /* Use local coordinate space */
    
    for (i=0; i<nbands; i++) Cx->coords[i] = tmp[i];

    Cx->nbands = nbands;
return 1;
}
/********************************************************************
                   outputting stats
********************************************************************/
/*********************************************
         DEF: print_CpMcnts()
Called after incremental build
*********************************************/
void print_CpMcnts()
{
  if (Cp.useFlag==0) { ZBuf[0] = '\0';
  }
  else {


  	/* Counts removed because it prevented parallelization */
    sprintf(ZBuf,
      "CpM counts: %d (%d %1.0e) %d (%d %1.0e) %d (%d %1.0e)\n", 
        -1, CpM[1].nm, CpM[1].cut,
        -1, CpM[2].nm, CpM[2].cut,
        -1, CpM[3].nm, CpM[3].cut);

  }
}

/*******************************************************
// test_write()
// Write out a file with all clone and matrix data in it
//
// sness
********************************************************/
void test_write(char *filename)
{
  int i, j;
  struct markertop *marker;
  FILE *outfile;
  FILE *hostfile;
  char line[256];
   
  hostfile = fopen("/etc/HOSTNAME","r");
  fgets(line,255,hostfile);
  strcat(filename,line);
  filename[strlen(filename)-1] = '\0';


  printf("test_write\n");fflush(stdout);
  if ((outfile = fopen(filename,"w")) == NULL) {
    printf("Could not open %s\n",filename); fflush(stdout);
 	exit(0);
  }

  /* Write out clones */
  fprintf(outfile,"Clonedata %i\n",arrayMax(acedata));
  fflush(outfile);
  for (i = 0; i < arrayMax(acedata); i++) {
 	fprintf(outfile,"Clone : %s\n",arr(acedata, i, CLONE).clone);
	fflush(outfile);

 	marker = arr(acedata, i, CLONE).marker;
 	j = 0;
 	while (marker) {
 	  fprintf(outfile,"Marker %i : %s\n",j,marker->marker);
	  fflush(outfile);
 	  marker = marker->nextmarker;
 	}
 	fprintf(outfile,"\n");
	fflush(outfile);

  }	

  /* Write out markers */
  fprintf(outfile,"Markerdata %i\n",arrayMax(markerdata));
  fflush(outfile);


  for (i = 0; i < arrayMax(markerdata); i++) {
 	fprintf(outfile,"Marker : %s\n",arr(markerdata, i, MARKER).marker);
	fflush(outfile);

  }

  /* Write out other information needed by ZscoreCpM */
  fprintf(outfile,"Pz.cutFlt = %.2e\n",Pz.cutFlt);
  fprintf(outfile,"Cp.useFlag = %i\n",Cp.useFlag);
  fprintf(outfile,"CpM[MAX_CpM-1].cut = %.2e\n",CpM[MAX_CpM-1].cut);
  fprintf(outfile,"MAX_CpM = %i\n",MAX_CpM);
  for (i = 0; i < MAX_CpM; i++) {
	fprintf(outfile,"CpM[%i].cut = %.2e\n",i,CpM[i].cut);
	fprintf(outfile,"CpM[%i].nm = %i\n",i,CpM[i].nm);
  }

  fflush(outfile);
  fclose(outfile);
}


/***********************************************************
                 DEF: fast_sulston
*************************************************************/
static void fast_sulston_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result)
{
register int i, k, n, nbH, nbL;
double a, c, pp, pa, pb;
static int first_time=1, tol = -1;
static double lim, Logfact[NBANDS+1], psmn=0.0;

    if (first_time) {
       Logfact[0] = 0.0;
       Logfact[1] = 0.0;
       for (i = 2; i <= NBANDS; ++i)
         Logfact[i] = Logfact[i - 1] + log((double) i);
       first_time=0;
       lim = log(DBL_MIN);
    }
    if (Pz.tol == 0)  psmn = 1 - (double) 1.0 / Proj.gel_len;
    else  psmn = 1 - (double) (Pz.tol << 1) / Proj.gel_len;
    tol = Pz.tol;
    
    nbL = MiN(clone1->nbands,clone2->nbands);
    nbH = MaX(clone1->nbands,clone2->nbands);

    result->prob = DBL_MIN;
    
    pp = (double) pow(psmn, nbH);
    pa = log(pp);
    pb = log(1.0 - pp);

    /* 18june00 starting at nbl 64, nhl 74, if the match is zero,
            this optimization causes the prob to be below cutoff */

    if (result->match==0) result->prob =1.0;
    else {
      for (k=0, n = result->match; n <= nbL && k < FAST_SUL_ITERATIONS; n++, k++)
      {
            c = Logfact[nbL] - Logfact[nbL - n] - Logfact[n];
            a = c + n * pb + (nbL - n) * pa;
            result->prob +=  (a >= lim) ? exp(a) : exp(lim);
       }
    }
}

/***********************************************************
                 DEF: pure_sulston
*************************************************************/
static void pure_sulston_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result)
{
register int i, k, n, nbH, nbL;
double a, c, pp, pa, pb;
static int first_time=1, tol=-1;
static double lim,  Logfact[NBANDS+1], psmn=0.0;

    if (first_time) {
       Logfact[0] = 0.0;
       Logfact[1] = 0.0;
       for (i = 2; i <= NBANDS; ++i)
         Logfact[i] = Logfact[i - 1] + log((double) i);
       first_time=0;
       lim = log(DBL_MIN);
    }
    if (tol!=Pz.tol) {
      if (Pz.tol == 0)  psmn = 1 - (double) 1.0 / Proj.gel_len;
      else  psmn = 1 - (double) (Pz.tol << 1) / Proj.gel_len;
      tol = Pz.tol;
    }
    nbL = MiN(clone1->nbands,clone2->nbands);
    nbH = MaX(clone1->nbands,clone2->nbands);

    result->prob = DBL_MIN;
    
    pp = (double) pow(psmn, nbH);
    pa = log(pp);
    pb = log(1.0 - pp);

    for (n = result->match; n <= nbL; n++)
    {
       c = Logfact[nbL] - Logfact[nbL - n] - Logfact[n];
       a = c + n * pb + (nbL - n) * pa;
       result->prob +=  (a >= lim) ? exp(a) : exp(lim);
    }
}


/***********************************************************
                 DEF: Zsulston
*************************************************************/
void Zsulston_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result)

{

    Zgetmatch_reentrant(clone1,clone2,result);

    if (0 == get_precomp_value(clone1->nbands, clone2->nbands, result->match, &result->prob))
    {
        if (Proj.eq2==1) Zmott_reentrant(clone1,clone2,result);
        else if (Proj.eq2==2 || clone1->nbands >=60 || clone2->nbands >= 60) 
             pure_sulston_reentrant(clone1,clone2,result);
        else fast_sulston_reentrant(clone1,clone2,result);
    }

} 

void choose_sulston(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result)
{
    if (Proj.eq2 == 1)
    {
       Zmott_reentrant(clone1,clone2,result);
    }
    else if (Proj.eq2==2 || clone1->nbands >=60 || clone2->nbands >= 60) 
    {
         pure_sulston_reentrant(clone1,clone2,result);
    }
    else 
    {
         fast_sulston_reentrant(clone1,clone2,result);
    }
}
/*********************************************************************
                    DEF: Zgetmatch

This accepts two clones and computes the number of matching bands
which it stores in result->

*******************************************************************/
void Zgetmatch_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result)
{
  int kbnd, k, l, idiff, lstart;
  
  result->total = result->match = result->newMatch = 0;
  lstart = 0;
  for (k = 0; k < clone1->nbands; ++k)
  {
      kbnd = clone1->coords[k];
      for (l = lstart; l < clone2->nbands; ++l)
      {
	  idiff = kbnd - clone2->coords[l];
	  if (Ztol(kbnd, idiff)) 
	    {
	      result->match++;
	      if (validDistribution) {
			/* trying different values for distribution */
	        result->newMatch += (distrib[kbnd] + distrib[kbnd-idiff])/2;
		result->match = result->newMatch;
	      }
	      lstart = l + 1;
	      break;
	    }
	  else if (idiff < 0)
	    {
	      lstart = l;
	      break;
	    }
       }
   }
  return;
} 

/**************************************************************
            DE: ZscoreCpM
All clam routines call this for scoreing (Zsulston or Zmott must be called
before this).
**************************************************************/
int ZscoreCpM_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result)
{
int junk, exponent=0;
char str[20];
register int i;

   result->mark = 0;
   if (result->prob < Pz.cutFlt) {
       CpMTbl[0].cnt++;
       sprintf(str,"%.0e", result->prob);
	   
       sscanf(str,"%de-%d",&junk, &exponent);

       if (exponent==0) {
		 printf("Clone1 %i\tScore%.0e\tScore%f\n", clone1->nbands,result->prob,result->prob);
		 
		 fflush(stdout);

		 printf("exp = 0\n");
	   }

       return exponent;
   }
   if (Cp.useFlag==0)  return 0;
   if (!Zcnt_shared_markers()) return 0;
   if (result->prob > CpM[MAX_CpM-1].cut) return 0;

   for (i=1; i<MAX_CpM; i++) 
      if (result->prob < CpM[i].cut && result->mark >= CpM[i].nm) {
         CpMTbl[i].cnt++;
         sprintf(str,"%.0e", result->prob);
         sscanf(str,"%de-%d",&junk, &exponent);
         if (exponent==0) printf("exp = 0\n");
         return exponent;
      }
   return 0;
}


/******************************************************************
                 DEF: Zcnt_shared_markers()
*****************************************************************/
int Zcnt_shared_markers_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result)
{
struct marker *markerptr;
struct markertop *mk1, *mk2;
CLONE *c1, *c2;
int score=0;

   c1 = arrp(acedata, clone1->cin, CLONE);
   c2 = arrp(acedata, clone2->cin, CLONE);
   for (mk1 = c1->marker; mk1 != NULL; mk1 = mk1->nextmarker)
   {
        markerptr = arrp(markerdata, mk1->markerindex, MARKER);
        if(!Cp.useYBP &&
             (markerptr->type == markYAC || markerptr->type == markBAC ||
                                markerptr->type == markPAC)) continue;
        if(!Cp.usePCR && markerptr->type == markPCR) continue;
        if(!Cp.useREP && markerptr->type == markREP) continue;

        for (mk2 = c2->marker; mk2 != NULL; mk2 = mk2->nextmarker)
            if (mk1->markerindex == mk2->markerindex) {
               score++;
               break;
            }
   }
   Sz.mark = score;
return score;
}

/***********************************************************
                 DEF: Zmott
***********************************************************/
void Zmott_reentrant(struct tmpCz* clone1, struct tmpCz* clone2, struct tmpSz* result)
{
double bin, b, lp1, lp2, p1, p2, p, tol, tol1, pt;

    p1 = (double) clone1->nbands / (double) Proj.gel_len;
    p2 = (double) clone2->nbands / (double) Proj.gel_len;
    p = p1*p2;
    lp1 = log(1.0-p1);
    lp2 = log(1.0-p2);
    tol = (double) Pz.tol;
    tol1 = (double) Pz.tol+1;
    pt = -p + p1*(1.0 - exp(tol1*lp2))*(exp(tol*lp2)) +
              p2*(1.0 - exp(tol1*lp1))*(exp(tol*lp1));
    if (validDistribution) {
      b = (double) result->newMatch;
    } else {
      b = (double) result->match;
    }

    bin = (double) Proj.gel_len;
    result->prob = (double) 
      exp(b*log(pt) + (bin-b)*log(1.0-pt) + 
      lgamma(bin+1.0) -lgamma(b+1.0) - lgamma(bin-b+1.0));
}

